///*
// * Licensed to the Apache Software Foundation (ASF) under one or more
// * contributor license agreements.  See the NOTICE file distributed with
// * this work for additional information regarding copyright ownership.
// * The ASF licenses this file to You under the Apache License, Version 2.0
// * (the "License"); you may not use this file except in compliance with
// * the License.  You may obtain a copy of the License at
// *
// *      http://www.apache.org/licenses/LICENSE-2.0
// *
// * Unless required by applicable law or agreed to in writing, software
// * distributed under the License is distributed on an "AS IS" BASIS,
// * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// * See the License for the specific language governing permissions and
// * limitations under the License.
// */
//
//package org.apache.commons.math4.legacy.linear;
//
//import org.apache.commons.numbers.complex.Complex;
//import org.apache.commons.numbers.core.Precision;
//import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
//import org.apache.commons.math4.legacy.exception.MathArithmeticException;
//import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
//import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
//import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
//import org.apache.commons.math4.core.jdkmath.JdkMath;
//
///**
// * Calculates the eigen decomposition of a real matrix.
// * <p>
// * The eigen decomposition of matrix A is a set of two matrices:
// * V and D such that A = V &times; D &times; V<sup>T</sup>.
// * A, V and D are all m &times; m matrices.
// * <p>
// * This class is similar in spirit to the {@code EigenvalueDecomposition}
// * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
// * library, with the following changes:
// * <ul>
// *   <li>a {@link #getVT() getVt} method has been added,</li>
// *   <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and
// *       {@link #getImagEigenvalue(int) getImagEigenvalue} methods to pick up a
// *       single eigenvalue have been added,</li>
// *   <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a
// *       single eigenvector has been added,</li>
// *   <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
// *   <li>a {@link #getSolver() getSolver} method has been added.</li>
// * </ul>
// * <p>
// * As of 3.1, this class supports general real matrices (both symmetric and non-symmetric):
// * <p>
// * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal
// * and the eigenvector matrix V is orthogonal, i.e.
// * {@code A = V.multiply(D.multiply(V.transpose()))} and
// * {@code V.multiply(V.transpose())} equals the identity matrix.
// * </p>
// * <p>
// * If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real
// * eigenvalues in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2
// * blocks:
// * <pre>
// *    [lambda, mu    ]
// *    [   -mu, lambda]
// * </pre>
// * The columns of V represent the eigenvectors in the sense that {@code A*V = V*D},
// * i.e. A.multiply(V) equals V.multiply(D).
// * The matrix V may be badly conditioned, or even singular, so the validity of the
// * equation {@code A = V*D*inverse(V)} depends upon the condition of V.
// * <p>
// * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
// * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
// * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
// * New-York.
// *
// * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
// * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
// * @since 2.0 (changed to concrete class in 3.0)
// */
//public class EigenDecomposition {
//    /** Internally used epsilon criteria. */
//    private static final double EPSILON = 1e-12;
//    /** Maximum number of iterations accepted in the implicit QL transformation. */
//    private static final byte MAX_ITER = 30;
//    /** Main diagonal of the tridiagonal matrix. */
//    private double[] main;
//    /** Secondary diagonal of the tridiagonal matrix. */
//    private double[] secondary;
//    /**
//     * Transformer to tridiagonal (may be null if matrix is already
//     * tridiagonal).
//     */
//    private TriDiagonalTransformer transformer;
//    /** Real part of the realEigenvalues. */
//    private double[] realEigenvalues;
//    /** Imaginary part of the realEigenvalues. */
//    private double[] imagEigenvalues;
//    /** Eigenvectors. */
//    private ArrayRealVector[] eigenvectors;
//    /** Cached value of V. */
//    private RealMatrix cachedV;
//    /** Cached value of D. */
//    private RealMatrix cachedD;
//    /** Cached value of Vt. */
//    private RealMatrix cachedVt;
//    /** Whether the matrix is symmetric. */
//    private final boolean isSymmetric;
//
//    /**
//     * Calculates the eigen decomposition of the given real matrix.
//     * <p>
//     * Supports decomposition of a general matrix since 3.1.
//     *
//     * @param matrix Matrix to decompose.
//     * @throws MaxCountExceededException if the algorithm fails to converge.
//     * @throws MathArithmeticException if the decomposition of a general matrix
//     * results in a matrix with zero norm
//     * @since 3.1
//     */
//    public EigenDecomposition(final RealMatrix matrix)
//        throws MathArithmeticException {
//        final double symTol = 10 * matrix.getRowDimension() * matrix.getColumnDimension() * Precision.EPSILON;
//        isSymmetric = MatrixUtils.isSymmetric(matrix, symTol);
//        if (isSymmetric) {
//            transformToTridiagonal(matrix);
//            findEigenVectors(transformer.getQ().getData());
//        } else {
//            final SchurTransformer t = transformToSchur(matrix);
//            findEigenVectorsFromSchur(t);
//        }
//    }
//
//    /**
//     * Calculates the eigen decomposition of the symmetric tridiagonal
//     * matrix.  The Householder matrix is assumed to be the identity matrix.
//     *
//     * @param main Main diagonal of the symmetric tridiagonal form.
//     * @param secondary Secondary of the tridiagonal form.
//     * @throws MaxCountExceededException if the algorithm fails to converge.
//     * @since 3.1
//     */
//    public EigenDecomposition(final double[] main, final double[] secondary) {
//        isSymmetric = true;
//        this.main      = main.clone();
//        this.secondary = secondary.clone();
//        transformer    = null;
//        final int size = main.length;
//        final double[][] z = new double[size][size];
//        for (int i = 0; i < size; i++) {
//            z[i][i] = 1.0;
//        }
//        findEigenVectors(z);
//    }
//
//    /**
//     * Gets the matrix V of the decomposition.
//     * V is an orthogonal matrix, i.e. its transpose is also its inverse.
//     * The columns of V are the eigenvectors of the original matrix.
//     * No assumption is made about the orientation of the system axes formed
//     * by the columns of V (e.g. in a 3-dimension space, V can form a left-
//     * or right-handed system).
//     *
//     * @return the V matrix.
//     */
//    public RealMatrix getV() {
//
//        if (cachedV == null) {
//            final int m = eigenvectors.length;
//            cachedV = MatrixUtils.createRealMatrix(m, m);
//            for (int k = 0; k < m; ++k) {
//                cachedV.setColumnVector(k, eigenvectors[k]);
//            }
//        }
//        // return the cached matrix
//        return cachedV;
//    }
//
//    /**
//     * Gets the block diagonal matrix D of the decomposition.
//     * D is a block diagonal matrix.
//     * Real eigenvalues are on the diagonal while complex values are on
//     * 2x2 blocks { {real +imaginary}, {-imaginary, real} }.
//     *
//     * @return the D matrix.
//     *
//     * @see #getRealEigenvalues()
//     * @see #getImagEigenvalues()
//     */
//    public RealMatrix getD() {
//
//        if (cachedD == null) {
//            // cache the matrix for subsequent calls
//            cachedD = MatrixUtils.createRealMatrixWithDiagonal(realEigenvalues);
//
//            for (int i = 0; i < imagEigenvalues.length; i++) {
//                if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) > 0) {
//                    cachedD.setEntry(i, i+1, imagEigenvalues[i]);
//                } else if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
//                    cachedD.setEntry(i, i-1, imagEigenvalues[i]);
//                }
//            }
//        }
//        return cachedD;
//    }
//
//    /**
//     * Gets the transpose of the matrix V of the decomposition.
//     * V is an orthogonal matrix, i.e. its transpose is also its inverse.
//     * The columns of V are the eigenvectors of the original matrix.
//     * No assumption is made about the orientation of the system axes formed
//     * by the columns of V (e.g. in a 3-dimension space, V can form a left-
//     * or right-handed system).
//     *
//     * @return the transpose of the V matrix.
//     */
//    public RealMatrix getVT() {
//
//        if (cachedVt == null) {
//            final int m = eigenvectors.length;
//            cachedVt = MatrixUtils.createRealMatrix(m, m);
//            for (int k = 0; k < m; ++k) {
//                cachedVt.setRowVector(k, eigenvectors[k]);
//            }
//        }
//
//        // return the cached matrix
//        return cachedVt;
//    }
//
//    /**
//     * Returns whether the calculated eigen values are complex or real.
//     * <p>The method performs a zero check for each element of the
//     * {@link #getImagEigenvalues()} array and returns {@code true} if any
//     * element is not equal to zero.
//     *
//     * @return {@code true} if the eigen values are complex, {@code false} otherwise
//     * @since 3.1
//     */
//    public boolean hasComplexEigenvalues() {
//        for (int i = 0; i < imagEigenvalues.length; i++) {
//            if (!Precision.equals(imagEigenvalues[i], 0.0, EPSILON)) {
//                return true;
//            }
//        }
//        return false;
//    }
//
//    /**
//     * Gets a copy of the real parts of the eigenvalues of the original matrix.
//     *
//     * @return a copy of the real parts of the eigenvalues of the original matrix.
//     *
//     * @see #getD()
//     * @see #getRealEigenvalue(int)
//     * @see #getImagEigenvalues()
//     */
//    public double[] getRealEigenvalues() {
//        return realEigenvalues.clone();
//    }
//
//    /**
//     * Returns the real part of the i<sup>th</sup> eigenvalue of the original
//     * matrix.
//     *
//     * @param i index of the eigenvalue (counting from 0)
//     * @return real part of the i<sup>th</sup> eigenvalue of the original
//     * matrix.
//     *
//     * @see #getD()
//     * @see #getRealEigenvalues()
//     * @see #getImagEigenvalue(int)
//     */
//    public double getRealEigenvalue(final int i) {
//        return realEigenvalues[i];
//    }
//
//    /**
//     * Gets a copy of the imaginary parts of the eigenvalues of the original
//     * matrix.
//     *
//     * @return a copy of the imaginary parts of the eigenvalues of the original
//     * matrix.
//     *
//     * @see #getD()
//     * @see #getImagEigenvalue(int)
//     * @see #getRealEigenvalues()
//     */
//    public double[] getImagEigenvalues() {
//        return imagEigenvalues.clone();
//    }
//
//    /**
//     * Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original
//     * matrix.
//     *
//     * @param i Index of the eigenvalue (counting from 0).
//     * @return the imaginary part of the i<sup>th</sup> eigenvalue of the original
//     * matrix.
//     *
//     * @see #getD()
//     * @see #getImagEigenvalues()
//     * @see #getRealEigenvalue(int)
//     */
//    public double getImagEigenvalue(final int i) {
//        return imagEigenvalues[i];
//    }
//
//    /**
//     * Gets a copy of the i<sup>th</sup> eigenvector of the original matrix.
//     *
//     * @param i Index of the eigenvector (counting from 0).
//     * @return a copy of the i<sup>th</sup> eigenvector of the original matrix.
//     * @see #getD()
//     */
//    public RealVector getEigenvector(final int i) {
//        return eigenvectors[i].copy();
//    }
//
//    /**
//     * Computes the determinant of the matrix.
//     *
//     * @return the determinant of the matrix.
//     */
//    public double getDeterminant() {
//        double determinant = 1;
//        for (double lambda : realEigenvalues) {
//            determinant *= lambda;
//        }
//        return determinant;
//    }
//
//    /**
//     * Computes the square-root of the matrix.
//     * This implementation assumes that the matrix is symmetric and positive
//     * definite.
//     *
//     * @return the square-root of the matrix.
//     * @throws MathUnsupportedOperationException if the matrix is not
//     * symmetric or not positive definite.
//     * @since 3.1
//     */
//    public RealMatrix getSquareRoot() {
//        if (!isSymmetric) {
//            throw new MathUnsupportedOperationException();
//        }
//
//        final double[] sqrtEigenValues = new double[realEigenvalues.length];
//        for (int i = 0; i < realEigenvalues.length; i++) {
//            final double eigen = realEigenvalues[i];
//            if (eigen <= 0) {
//                throw new MathUnsupportedOperationException();
//            }
//            sqrtEigenValues[i] = JdkMath.sqrt(eigen);
//        }
//        final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues);
//        final RealMatrix v = getV();
//        final RealMatrix vT = getVT();
//
//        return v.multiply(sqrtEigen).multiply(vT);
//    }
//
//    /**
//     * Gets a solver for finding the A &times; X = B solution in exact
//     * linear sense.
//     * <p>
//     * Since 3.1, eigen decomposition of a general matrix is supported,
//     * but the {@link DecompositionSolver} only supports real eigenvalues.
//     *
//     * @return a solver
//     * @throws MathUnsupportedOperationException if the decomposition resulted in
//     * complex eigenvalues
//     */
//    public DecompositionSolver getSolver() {
//        if (hasComplexEigenvalues()) {
//            throw new MathUnsupportedOperationException();
//        }
//        return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
//    }
//
//    /** Specialized solver. */
//    private static final class Solver implements DecompositionSolver {
//        /** Real part of the realEigenvalues. */
//        private final double[] realEigenvalues;
//        /** Imaginary part of the realEigenvalues. */
//        private final double[] imagEigenvalues;
//        /** Eigenvectors. */
//        private final ArrayRealVector[] eigenvectors;
//
//        /**
//         * Builds a solver from decomposed matrix.
//         *
//         * @param realEigenvalues Real parts of the eigenvalues.
//         * @param imagEigenvalues Imaginary parts of the eigenvalues.
//         * @param eigenvectors Eigenvectors.
//         */
//        private Solver(final double[] realEigenvalues,
//                final double[] imagEigenvalues,
//                final ArrayRealVector[] eigenvectors) {
//            this.realEigenvalues = realEigenvalues;
//            this.imagEigenvalues = imagEigenvalues;
//            this.eigenvectors = eigenvectors;
//        }
//
//        /**
//         * Solves the linear equation A &times; X = B for symmetric matrices A.
//         * <p>
//         * This method only finds exact linear solutions, i.e. solutions for
//         * which ||A &times; X - B|| is exactly 0.
//         * </p>
//         *
//         * @param b Right-hand side of the equation A &times; X = B.
//         * @return a Vector X that minimizes the two norm of A &times; X - B.
//         *
//         * @throws DimensionMismatchException if the matrices dimensions do not match.
//         * @throws SingularMatrixException if the decomposed matrix is singular.
//         */
//        @Override
//        public RealVector solve(final RealVector b) {
//            if (!isNonSingular()) {
//                throw new SingularMatrixException();
//            }
//
//            final int m = realEigenvalues.length;
//            if (b.getDimension() != m) {
//                throw new DimensionMismatchException(b.getDimension(), m);
//            }
//
//            final double[] bp = new double[m];
//            for (int i = 0; i < m; ++i) {
//                final ArrayRealVector v = eigenvectors[i];
//                final double[] vData = v.getDataRef();
//                final double s = v.dotProduct(b) / realEigenvalues[i];
//                for (int j = 0; j < m; ++j) {
//                    bp[j] += s * vData[j];
//                }
//            }
//
//            return new ArrayRealVector(bp, false);
//        }
//
//        /** {@inheritDoc} */
//        @Override
//        public RealMatrix solve(RealMatrix b) {
//
//            if (!isNonSingular()) {
//                throw new SingularMatrixException();
//            }
//
//            final int m = realEigenvalues.length;
//            if (b.getRowDimension() != m) {
//                throw new DimensionMismatchException(b.getRowDimension(), m);
//            }
//
//            final int nColB = b.getColumnDimension();
//            final double[][] bp = new double[m][nColB];
//            final double[] tmpCol = new double[m];
//            for (int k = 0; k < nColB; ++k) {
//                for (int i = 0; i < m; ++i) {
//                    tmpCol[i] = b.getEntry(i, k);
//                    bp[i][k]  = 0;
//                }
//                for (int i = 0; i < m; ++i) {
//                    final ArrayRealVector v = eigenvectors[i];
//                    final double[] vData = v.getDataRef();
//                    double s = 0;
//                    for (int j = 0; j < m; ++j) {
//                        s += v.getEntry(j) * tmpCol[j];
//                    }
//                    s /= realEigenvalues[i];
//                    for (int j = 0; j < m; ++j) {
//                        bp[j][k] += s * vData[j];
//                    }
//                }
//            }
//
//            return new Array2DRowRealMatrix(bp, false);
//        }
//
//        /**
//         * Checks whether the decomposed matrix is non-singular.
//         *
//         * @return true if the decomposed matrix is non-singular.
//         */
//        @Override
//        public boolean isNonSingular() {
//            double largestEigenvalueNorm = 0.0;
//            // Looping over all values (in case they are not sorted in decreasing
//            // order of their norm).
//            for (int i = 0; i < realEigenvalues.length; ++i) {
//                largestEigenvalueNorm = JdkMath.max(largestEigenvalueNorm, eigenvalueNorm(i));
//            }
//            // Corner case: zero matrix, all exactly 0 eigenvalues
//            if (largestEigenvalueNorm == 0.0) {
//                return false;
//            }
//            for (int i = 0; i < realEigenvalues.length; ++i) {
//                // Looking for eigenvalues that are 0, where we consider anything much much smaller
//                // than the largest eigenvalue to be effectively 0.
//                if (Precision.equals(eigenvalueNorm(i) / largestEigenvalueNorm, 0, EPSILON)) {
//                    return false;
//                }
//            }
//            return true;
//        }
//
//        /**
//         * @param i which eigenvalue to find the norm of
//         * @return the norm of ith (complex) eigenvalue.
//         */
//        private double eigenvalueNorm(int i) {
//            final double re = realEigenvalues[i];
//            final double im = imagEigenvalues[i];
//            return JdkMath.sqrt(re * re + im * im);
//        }
//
//        /**
//         * Get the inverse of the decomposed matrix.
//         *
//         * @return the inverse matrix.
//         * @throws SingularMatrixException if the decomposed matrix is singular.
//         */
//        @Override
//        public RealMatrix getInverse() {
//            if (!isNonSingular()) {
//                throw new SingularMatrixException();
//            }
//
//            final int m = realEigenvalues.length;
//            final double[][] invData = new double[m][m];
//
//            for (int i = 0; i < m; ++i) {
//                final double[] invI = invData[i];
//                for (int j = 0; j < m; ++j) {
//                    double invIJ = 0;
//                    for (int k = 0; k < m; ++k) {
//                        final double[] vK = eigenvectors[k].getDataRef();
//                        invIJ += vK[i] * vK[j] / realEigenvalues[k];
//                    }
//                    invI[j] = invIJ;
//                }
//            }
//            return MatrixUtils.createRealMatrix(invData);
//        }
//    }
//
//    /**
//     * Transforms the matrix to tridiagonal form.
//     *
//     * @param matrix Matrix to transform.
//     */
//    private void transformToTridiagonal(final RealMatrix matrix) {
//        // transform the matrix to tridiagonal
//        transformer = new TriDiagonalTransformer(matrix);
//        main = transformer.getMainDiagonalRef();
//        secondary = transformer.getSecondaryDiagonalRef();
//    }
//
//    /**
//     * Find eigenvalues and eigenvectors (Dubrulle et al., 1971).
//     *
//     * @param householderMatrix Householder matrix of the transformation
//     * to tridiagonal form.
//     */
//    private void findEigenVectors(final double[][] householderMatrix) {
//        final double[][]z = householderMatrix.clone();
//        final int n = main.length;
//        realEigenvalues = new double[n];
//        imagEigenvalues = new double[n];
//        final double[] e = new double[n];
//        for (int i = 0; i < n - 1; i++) {
//            realEigenvalues[i] = main[i];
//            e[i] = secondary[i];
//        }
//        realEigenvalues[n - 1] = main[n - 1];
//        e[n - 1] = 0;
//
//        // Determine the largest main and secondary value in absolute term.
//        double maxAbsoluteValue = 0;
//        for (int i = 0; i < n; i++) {
//            if (JdkMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
//                maxAbsoluteValue = JdkMath.abs(realEigenvalues[i]);
//            }
//            if (JdkMath.abs(e[i]) > maxAbsoluteValue) {
//                maxAbsoluteValue = JdkMath.abs(e[i]);
//            }
//        }
//        // Make null any main and secondary value too small to be significant
//        if (maxAbsoluteValue != 0) {
//            for (int i=0; i < n; i++) {
//                if (JdkMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
//                    realEigenvalues[i] = 0;
//                }
//                if (JdkMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
//                    e[i]=0;
//                }
//            }
//        }
//
//        for (int j = 0; j < n; j++) {
//            int its = 0;
//            int m;
//            do {
//                for (m = j; m < n - 1; m++) {
//                    double delta = JdkMath.abs(realEigenvalues[m]) +
//                        JdkMath.abs(realEigenvalues[m + 1]);
//                    if (JdkMath.abs(e[m]) + delta == delta) {
//                        break;
//                    }
//                }
//                if (m != j) {
//                    if (its == MAX_ITER) {
//                        throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
//                                                            MAX_ITER);
//                    }
//                    its++;
//                    double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
//                    double t = JdkMath.sqrt(1 + q * q);
//                    if (q < 0.0) {
//                        q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
//                    } else {
//                        q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
//                    }
//                    double u = 0.0;
//                    double s = 1.0;
//                    double c = 1.0;
//                    int i;
//                    for (i = m - 1; i >= j; i--) {
//                        double p = s * e[i];
//                        double h = c * e[i];
//                        if (JdkMath.abs(p) >= JdkMath.abs(q)) {
//                            c = q / p;
//                            t = JdkMath.sqrt(c * c + 1.0);
//                            e[i + 1] = p * t;
//                            s = 1.0 / t;
//                            c *= s;
//                        } else {
//                            s = p / q;
//                            t = JdkMath.sqrt(s * s + 1.0);
//                            e[i + 1] = q * t;
//                            c = 1.0 / t;
//                            s *= c;
//                        }
//                        if (e[i + 1] == 0.0) {
//                            realEigenvalues[i + 1] -= u;
//                            e[m] = 0.0;
//                            break;
//                        }
//                        q = realEigenvalues[i + 1] - u;
//                        t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
//                        u = s * t;
//                        realEigenvalues[i + 1] = q + u;
//                        q = c * t - h;
//                        for (int ia = 0; ia < n; ia++) {
//                            p = z[ia][i + 1];
//                            z[ia][i + 1] = s * z[ia][i] + c * p;
//                            z[ia][i] = c * z[ia][i] - s * p;
//                        }
//                    }
//                    if (t == 0.0 && i >= j) {
//                        continue;
//                    }
//                    realEigenvalues[j] -= u;
//                    e[j] = q;
//                    e[m] = 0.0;
//                }
//            } while (m != j);
//        }
//
//        //Sort the eigen values (and vectors) in increase order
//        for (int i = 0; i < n; i++) {
//            int k = i;
//            double p = realEigenvalues[i];
//            for (int j = i + 1; j < n; j++) {
//                if (realEigenvalues[j] > p) {
//                    k = j;
//                    p = realEigenvalues[j];
//                }
//            }
//            if (k != i) {
//                realEigenvalues[k] = realEigenvalues[i];
//                realEigenvalues[i] = p;
//                for (int j = 0; j < n; j++) {
//                    p = z[j][i];
//                    z[j][i] = z[j][k];
//                    z[j][k] = p;
//                }
//            }
//        }
//
//        // Determine the largest eigen value in absolute term.
//        maxAbsoluteValue = 0;
//        for (int i = 0; i < n; i++) {
//            if (JdkMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
//                maxAbsoluteValue=JdkMath.abs(realEigenvalues[i]);
//            }
//        }
//        // Make null any eigen value too small to be significant
//        if (maxAbsoluteValue != 0.0) {
//            for (int i=0; i < n; i++) {
//                if (JdkMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
//                    realEigenvalues[i] = 0;
//                }
//            }
//        }
//        eigenvectors = new ArrayRealVector[n];
//        final double[] tmp = new double[n];
//        for (int i = 0; i < n; i++) {
//            for (int j = 0; j < n; j++) {
//                tmp[j] = z[j][i];
//            }
//            eigenvectors[i] = new ArrayRealVector(tmp);
//        }
//    }
//
//    /**
//     * Transforms the matrix to Schur form and calculates the eigenvalues.
//     *
//     * @param matrix Matrix to transform.
//     * @return the {@link SchurTransformer Shur transform} for this matrix
//     */
//    private SchurTransformer transformToSchur(final RealMatrix matrix) {
//        final SchurTransformer schurTransform = new SchurTransformer(matrix);
//        final double[][] matT = schurTransform.getT().getData();
//
//        realEigenvalues = new double[matT.length];
//        imagEigenvalues = new double[matT.length];
//
//        for (int i = 0; i < realEigenvalues.length; i++) {
//            if (i == (realEigenvalues.length - 1) ||
//                Precision.equals(matT[i + 1][i], 0.0, EPSILON)) {
//                realEigenvalues[i] = matT[i][i];
//            } else {
//                final double x = matT[i + 1][i + 1];
//                final double p = 0.5 * (matT[i][i] - x);
//                final double z = JdkMath.sqrt(JdkMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1]));
//                realEigenvalues[i] = x + p;
//                imagEigenvalues[i] = z;
//                realEigenvalues[i + 1] = x + p;
//                imagEigenvalues[i + 1] = -z;
//                i++;
//            }
//        }
//        return schurTransform;
//    }
//
//    /**
//     * Performs a division of two complex numbers.
//     *
//     * @param xr real part of the first number
//     * @param xi imaginary part of the first number
//     * @param yr real part of the second number
//     * @param yi imaginary part of the second number
//     * @return result of the complex division
//     */
//    private Complex cdiv(final double xr, final double xi,
//                         final double yr, final double yi) {
//        return Complex.ofCartesian(xr, xi).divide(Complex.ofCartesian(yr, yi));
//    }
//
//    /**
//     * Find eigenvectors from a matrix transformed to Schur form.
//     *
//     * @param schur the schur transformation of the matrix
//     * @throws MathArithmeticException if the Schur form has a norm of zero
//     */
//    private void findEigenVectorsFromSchur(final SchurTransformer schur)
//        throws MathArithmeticException {
//        final double[][] matrixT = schur.getT().getData();
//        final double[][] matrixP = schur.getP().getData();
//
//        final int n = matrixT.length;
//
//        // compute matrix norm
//        double norm = 0.0;
//        for (int i = 0; i < n; i++) {
//           for (int j = JdkMath.max(i - 1, 0); j < n; j++) {
//               norm += JdkMath.abs(matrixT[i][j]);
//           }
//        }
//
//        // we can not handle a matrix with zero norm
//        if (Precision.equals(norm, 0.0, EPSILON)) {
//           throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
//        }
//
//        // Backsubstitute to find vectors of upper triangular form
//
//        double r = 0.0;
//        double s = 0.0;
//        double z = 0.0;
//
//        for (int idx = n - 1; idx >= 0; idx--) {
//            double p = realEigenvalues[idx];
//            double q = imagEigenvalues[idx];
//
//            if (Precision.equals(q, 0.0)) {
//                // Real vector
//                int l = idx;
//                matrixT[idx][idx] = 1.0;
//                for (int i = idx - 1; i >= 0; i--) {
//                    double w = matrixT[i][i] - p;
//                    r = 0.0;
//                    for (int j = l; j <= idx; j++) {
//                        r += matrixT[i][j] * matrixT[j][idx];
//                    }
//                    if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
//                        z = w;
//                        s = r;
//                    } else {
//                        l = i;
//                        if (Precision.equals(imagEigenvalues[i], 0.0)) {
//                            if (w != 0.0) {
//                                matrixT[i][idx] = -r / w;
//                            } else {
//                                matrixT[i][idx] = -r / (Precision.EPSILON * norm);
//                            }
//                        } else {
//                            // Solve real equations
//                            double x = matrixT[i][i + 1];
//                            double y = matrixT[i + 1][i];
//                            q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
//                                imagEigenvalues[i] * imagEigenvalues[i];
//                            double t = (x * s - z * r) / q;
//                            matrixT[i][idx] = t;
//                            if (JdkMath.abs(x) > JdkMath.abs(z)) {
//                                matrixT[i + 1][idx] = (-r - w * t) / x;
//                            } else {
//                                matrixT[i + 1][idx] = (-s - y * t) / z;
//                            }
//                        }
//
//                        // Overflow control
//                        double t = JdkMath.abs(matrixT[i][idx]);
//                        if ((Precision.EPSILON * t) * t > 1) {
//                            for (int j = i; j <= idx; j++) {
//                                matrixT[j][idx] /= t;
//                            }
//                        }
//                    }
//                }
//            } else if (q < 0.0) {
//                // Complex vector
//                int l = idx - 1;
//
//                // Last vector component imaginary so matrix is triangular
//                if (JdkMath.abs(matrixT[idx][idx - 1]) > JdkMath.abs(matrixT[idx - 1][idx])) {
//                    matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
//                    matrixT[idx - 1][idx]     = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
//                } else {
//                    final Complex result = cdiv(0.0, -matrixT[idx - 1][idx],
//                                                matrixT[idx - 1][idx - 1] - p, q);
//                    matrixT[idx - 1][idx - 1] = result.getReal();
//                    matrixT[idx - 1][idx]     = result.getImaginary();
//                }
//
//                matrixT[idx][idx - 1] = 0.0;
//                matrixT[idx][idx]     = 1.0;
//
//                for (int i = idx - 2; i >= 0; i--) {
//                    double ra = 0.0;
//                    double sa = 0.0;
//                    for (int j = l; j <= idx; j++) {
//                        ra += matrixT[i][j] * matrixT[j][idx - 1];
//                        sa += matrixT[i][j] * matrixT[j][idx];
//                    }
//                    double w = matrixT[i][i] - p;
//
//                    if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
//                        z = w;
//                        r = ra;
//                        s = sa;
//                    } else {
//                        l = i;
//                        if (Precision.equals(imagEigenvalues[i], 0.0)) {
//                            final Complex c = cdiv(-ra, -sa, w, q);
//                            matrixT[i][idx - 1] = c.getReal();
//                            matrixT[i][idx] = c.getImaginary();
//                        } else {
//                            // Solve complex equations
//                            double x = matrixT[i][i + 1];
//                            double y = matrixT[i + 1][i];
//                            double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
//                                        imagEigenvalues[i] * imagEigenvalues[i] - q * q;
//                            final double vi = (realEigenvalues[i] - p) * 2.0 * q;
//                            if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
//                                vr = Precision.EPSILON * norm *
//                                     (JdkMath.abs(w) + JdkMath.abs(q) + JdkMath.abs(x) +
//                                      JdkMath.abs(y) + JdkMath.abs(z));
//                            }
//                            final Complex c     = cdiv(x * r - z * ra + q * sa,
//                                                       x * s - z * sa - q * ra, vr, vi);
//                            matrixT[i][idx - 1] = c.getReal();
//                            matrixT[i][idx]     = c.getImaginary();
//
//                            if (JdkMath.abs(x) > (JdkMath.abs(z) + JdkMath.abs(q))) {
//                                matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] +
//                                                           q * matrixT[i][idx]) / x;
//                                matrixT[i + 1][idx]     = (-sa - w * matrixT[i][idx] -
//                                                           q * matrixT[i][idx - 1]) / x;
//                            } else {
//                                final Complex c2        = cdiv(-r - y * matrixT[i][idx - 1],
//                                                               -s - y * matrixT[i][idx], z, q);
//                                matrixT[i + 1][idx - 1] = c2.getReal();
//                                matrixT[i + 1][idx]     = c2.getImaginary();
//                            }
//                        }
//
//                        // Overflow control
//                        double t = JdkMath.max(JdkMath.abs(matrixT[i][idx - 1]),
//                                                JdkMath.abs(matrixT[i][idx]));
//                        if ((Precision.EPSILON * t) * t > 1) {
//                            for (int j = i; j <= idx; j++) {
//                                matrixT[j][idx - 1] /= t;
//                                matrixT[j][idx] /= t;
//                            }
//                        }
//                    }
//                }
//            }
//        }
//
//        // Back transformation to get eigenvectors of original matrix
//        for (int j = n - 1; j >= 0; j--) {
//            for (int i = 0; i <= n - 1; i++) {
//                z = 0.0;
//                for (int k = 0; k <= JdkMath.min(j, n - 1); k++) {
//                    z += matrixP[i][k] * matrixT[k][j];
//                }
//                matrixP[i][j] = z;
//            }
//        }
//
//        eigenvectors = new ArrayRealVector[n];
//        final double[] tmp = new double[n];
//        for (int i = 0; i < n; i++) {
//            for (int j = 0; j < n; j++) {
//                tmp[j] = matrixP[j][i];
//            }
//            eigenvectors[i] = new ArrayRealVector(tmp);
//        }
//    }
//}
